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Abstract 

We test a recent proposal to use approximate trivializing maps in a 
field theory to speed up Hybrid Monte Carlo simulations. Simulating 
the CP^ -1 model, we find a small improvement with the leading order 
transformation, which is however compensated by the additional compu- 
tational overhead. The scaling of the algorithm towards the continuum 
is not changed. In particular, the effect of the topological modes on the 
autocorrelation times is studied. 

1 Introduction 

In the simulation of statistical models, many Monte Carlo methods experience a 
significant increase in effort when approaching a continuous phase transition of 
the theory. This phenomenon is called critical slowing down and depends strongly 
on the nature of the underlying theory and the algorithm used. For some models, 
algorithms have been found which completely eliminate this slowing down or even 
lead to a speed-up when the critical line is approached. A particular type of 
critical slowing down is associated with the topological modes of the theories. 
In QCD, e.g., the topological charge of the gauge configuration is known to be 
particularly problematic with both single link updates [1] and algorithms based 
on molecular dynamics [2]. 

Recently, a possible solution to the problem has been proposed [3] , for which 
field transformations given by flow equations are introduced. Exactly integrat- 
ing the flow equations, the theory becomes trivial and therefore also trivial to 
simulate. The problem is that the differential equations generating the flow are 
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not exactly known, however, the first terms of a power series of the correspond- 
ing kernel can easily be constructed. At this point, it is unclear whether it is 
sufficient to just know the first order of the differential equations to sufficiently 
mitigate the problem. Also the accuracy to which the differential equations have 
to be integrated is unknown. 

In this letter, we therefore put this proposal to a test. Since it should work 
for any field theory, we simulate the CP 7V_1 model for N = 10 using the Hybrid 
Monte Carlo (HMC) algorithm [3] , which is an integral part of the program laid 
out in Ref. [3J. We choose this model, because it shares some similarities with 
QCD like asymptotic freedom and confinement and, like in QCD, the topological 
modes show a much more severe critical slowing down than other observables [5] . 
Most importantly, the accuracy achievable in this two dimensional model is much 
higher than in QCD and the reduced cost also allows for a better mapping of the 
rather high dimensional parameter space of the problem. 

In Sect. [2] we give the details of the model and how to set up the HMC 
algorithm for it. After that, we discuss the trivializing map and its construction 
to leading order in Sec. |3j Then we give the parameters of the simulation which 
enables us to put the approximate trivializing map to a test whose results we 
give in Sect. |4| In the definition of the action, the observables and as a point of 
reference for the main quantities, we follow the paper by Campostrini, Rossi and 
Vicari [6]. 

2 The CP^model and Hybrid Monte Carlo 

We immediately give the model on a square lattice with lattice spacing a and 
sites n = a(ni,ri2), U\ and 712 integer numbers. The fields living on the sites 
are complex N component unit vectors z n connected by U(l) link variables A n>At , 
which are represented by complex numbers on the unit circle. The action is then 
given by [7] 

2 

s[z, \] = -npJ2J2 ( z l+fi z nK, + 4 - 2) . (i) 

n fJ.=l 

The gauge fields X n>fl can be integrated out analytically, thus this lattice action 
is expected to lie within the universality class of the CP Ar_1 model. 

For the Hybrid Monte Carlo algorithm as well as for the field transformation 
we will need derivatives with respect to the field degrees of freedom. Therefore 
it is convenient to treat the fields as real 2N component fields x n , which live on 
the unit sphere in M. 2N 

X n ,2i = Re(>n,i) , Xn,2i+1 = W^j) , % = 0, . . . , JV - 1 . (2) 

In the following, we will make use of either z or x, such that each formula appears 
in its most simple form. 
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The U(l) fields A are naturally parametrized by the angle G [0, 2ir) with 
A = e tcf> . Its action on the vectors x is then given by SO(2) 2N x 2N matrices 
A n , M (0), which are zero everywhere but on the diagonal 2x2 blocks 

( A 2 i,2i A 2i ,2i+i \ = /cos -sin (A . = Q N _ 1 
\Mi+i,2i A 2 i+i,2i+i/ \sm(f) cos0 J ' 

With these definitions, we can write the action as function of real variables only 

S\x^]=-N^{ xT n J n-^) , (3) 

n 

where we introduced the "spin sum" of gauge-transported nearest neighbors 

±2 

J n = A L X «+A > with A «,-m : = A I-a,m • ( 4 ) 

|U=±1 

The partition function can then be easily written by embedding the unit vectors 
x n into R 2N . 

2.1 Observables 

We will focus our study on a few, central observables. Mainly the energy density 
E = S/(N/3V), the magnetic susceptibility \m and the correlation length £q. 
The latter two are constructed from the two point function in momentum space 

G P (k) = ^ y^(trP n P m ) conn exp f-p(" -m)-kj (5) 
with P n = ZnZ^. Using these definitions, the two remaining observables are then 



The topological charge density g„ is given by the sum over the angles between 
the spins around a plaquette 

Qn = -^tnv{9n,n + 0n+fi,u ~ 6 n +u,u ~ Q n ,v) m od 1 ; — - < q n < - 

with 9 n ^ = axg(zj l z n+ jj) . The topological charge is the volume integral of this 
quantity <? = En Qn- 
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2.2 Hybrid Monte Carlo 



In Hybrid Monte Carlo, the fields x and are updated by introducing conjugate 
momenta it and u and solving classical equations of motion associated with the 
Hamiltonian 

H[ir, u, x, 4>] = - ^2(ir n ) 2 + - ^(w„ iM ) 2 + S[x, 4>] . 

n n,fi 

The momenta live in the tangent spaces of the respective field manifolds and 
therefore ir n G M. 2N with 7r n • x n = and u n ^ G K without further conditions, 
because the manifold is flat. The Hamilton equations of motion are then 

n n = -V Xn S[x, (j)} , = -d^S[x, <p\ . (8) 

A natural derivative d l x of a function f(x) defined on the unit sphere is the 
projection of the ordinary gradient V x in M? N onto the tangent space of the 
sphere 

dif(x)=[(l-xx T )V x f(x)].. (9) 

This corresponds to continuing the function f(x) to the full M. 2N via f(x) = 
f(x/\x\) and then taking ordinary derivatives. The forces in the equations of 
motion ^ then read for the action given in Eq. ^ 

K = ~V Xn S[x, 0] = 2N/3 (1 - x n x T n ) J n =: 2Nf3p n , (10) 

where we have defined p n as the projection of the spin sum J n to the tangent 
space at x n . The forces F$ for the conjugate momenta uj n ^ are 

F^ = -d^S[x^] = -2NPx T n TA^x n+fl with r = (11) 

The T is the translation of the imaginary unit i to the language of the 2 component 
real vectors. 

In the numerical simulations, we use a leap-frog integration scheme with a 
single time scale. For this, we need finite step size updates of the fields to numer- 
ically solve Eqs. (J8|. The only non trivial part is the update of the field x n with 
momentum 7r n , because also the updated variables have to fulfill the constraints 
\x' n \ = 1 and x' n ■ 7c' n = 0. For an infinitesimal step of size e, we therefore use the 
map $ e 

( X ') = = ( n a ^ Sina ) h with a = e\n\. (12) 

\7r J y— |7r|sma cos a / V 71 "/ 

It corresponds to the exact solution of the equation of motion in the absence of 
the forces F but subject to the constraint \x\ = 1. 
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3 Trivializing map in the CP^ ^odel 

The goal of the field trivialization is to find a map J 7 in field space such that 
the Jacobian J of the transformation compensates the action. For the partition 
function Eq. ([T| this would mean a transformation (x' } 0') = J r ^ 1 (x, 0) such that 

[dx][dcj>]e- s ^ = [ [ dx ^[ d 4>] e -sW,#)]+i oe a*j[s>j] (13) 



with the exponent equal to a constant. Since in this case all configurations in the 
new variables are equally likely, the molecular dynamics evolution of the HMC 
algorithm would not experience any forces and be very efficient. But also in the 
situation that one can only find an approximation to the exact trivializing map 
J 7 , one can expect a significant gain in the performance of the algorithm. 

In our simulations, the forces associated to the U(l) field A are much smaller 
than those associated to the spin variables z. For all considered — 0.7, . . . , 1.0 
in CP 9 , we found the ratio of the average forces to be about ten and also for the 
maximal force a factor of almost four. We therefore perform the field transfor- 
mation only on the x, leaving the <fi untouched. In the remaining part of this 
section, we go through the major steps of the computation, following the lines of 
Ref. [3]. 

The trivializing map J 7 can be obtained by integrating a flow T from t — to 
t — 1. Note, however, that integration to some tx < 1 will probably be a better 
choice for only approximate trivializing flows. A possible ansatz for the flow T is 
to take a gradient of an action S 

x* n (t) = -dl l S[t,x(t)]=Ti l [t,x(t)} . (14) 

The action S can be expanded in a power series in t 



S[t,x] = J2 tk S (k) l 



x\ 

k=0 



for which the leading order term can be constructed easily. One result of Ref. [3] 
is that the leading term of 5* fulfills 



Using the derivative defined in Eq. (|9j) , it is easy to show that on the unit sphere 

-&&f{x) = (2N - l)x ■ Vf(x) - tr [(1 - xx T )H f (x)] 
with the Hessian Hf(x)ij = didjf(x). This immediately leads to 

§m = 2W^T) S ■ (15) 
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The corresponding leading order trivializing flow from Eq. ( 14 ) reads 



T(°) = ™l-p (16) 
2N-l Pn ' 1 j 



As next step, we need a numerical integration scheme for the flow equation (14). 



Following [3] , we use an Euler integrator in which each step is similar to the finite 



step-size update discussed in Sec. 2.2 



TJt) 

x n {t + e s ) = cosa n x(t) + sina n " with a n = e s \T n \ . (17) 

This integrator has 0(e s ) errors when integrating to a fixed t, but as we see 
below, the method does not suffer significantly from these inaccuracies. To get the 
action in the transformed variables, the determinant of the Jacobian of the field 



transformation has to be computed, see Eq. (13). Since this is too complicated 
if all spins are changed at once, we follow Ref. [3] again and transform one spin 
at a time, sweeping through the lattice. This sweep is done for each step in e s of 
the Euler integrator. 



For the transformation given in Eq. (17), the determinant can be easily com- 



puted. In the language of Eq. (13), for a single step the primed quantities are at 
t, whereas the unprimed quantities at t + e s . This transformation only changes 
the angle 9 between the neighbor sum J n and the transformed spin x n , leaving 
the components perpendicular to this plane untouched. Specifically, it is changed 
to 

9 — 9' — oi — 9' — 1 J n \ sin9' . 

2N — 1 1 

Since the integration measure for the angular component of spherical coordinates 
in the M. 2N is (sin6 l ) 2Ar_3 (icos6 l one easily obtains the Jacobi determinant as 



det J n = ( 1 — - * Jn x 'n) ( COS a V~~T ^ n a 



2N-1 n n J \ \p' r 

The forces corresponding to the action constructed from the smoothed fields x(t), 
have to be computed using the chain rule. Since this is a standard procedure, we 
do not describe it here in detail. 



4 Details of the simulation 

To our knowledge, this investigation is the first using the HMC algorithm to 
simulate the CP Ar_1 model. It is certainly not the algorithm of choice for this 
theory, but our objective is the study of the improvement in the algorithm brought 
by the leading order trivialization. For comparison with the literature, we relied 
heavily on Ref. [6j, from which we reproduced several observables for all parameter 
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sets we looked at. Since there is no prior experience with the HMC in this model, 
we have to first study it in the normal variables and can then assess the change 
that is brought by the field transformation. We use the acronym THMC for the 
HMC with field transformation in the following. 

The figure of merit is the autocorrelation time of interesting observables in 
units of the molecular dynamics time. It is defined via the autocorrelation func- 
tion 

T A (t) = ((A(s + t)-A)(A(s)-A)) , 

where the average is over independent realizations of the Markov chain and A is 
the expectation value of the observable A. The integrated autocorrelation time 
is then 



Tint{A) = 1 - + y • (19) 

mtv ) 2 ^r A (o) v ; 



Since Monte Carlo histories are never infinitely long, the sum in Eq. ( 19 ) has to 
be truncated at some window W. For its choice, one has to balance the statis- 
tical error, which increases with larger windows, with the systematic error from 
neglecting the contributions beyond W . For this purpose, we use the software 
described in Ref. [8]. In most cases, its automatic criterion turned out to be 
sufficient, however, due to our high statistics data, we sometimes had to increase 
the parameter S, which influences the relative size between the window W cho- 
sen and ri n t(W). Also because of the high statistics, we are confident that the 
systematic error due to slow modes is under control in all our data points. In 
the following, all autocorrelation times are given in units of molecular dynamics 
time. 



4.1 Tuning of the parameters 

The HMC algorithm has two tuning parameters, the trajectory length r tra j and 
the accuracy with which the molecular dynamics equations are integrated. With 
the field transformation, one also has to fix the integration length tx of the field 
transformation and its step size. Let us go through these parameters, the final 
values are listed in Table Q3 

4.1.1 Trajectory length 

The effect of the trajectory length on the autocorrelation times for the CP 9 
model with /3 = 0.7 can be found in Fig. [TJ The left hand plot shows the HMC 
without trivialization, the right hand side the THMC with the leading order 
flow integrated up to tx = 0.47 (such that the force is minimized, as will be 
discussed later) with one Euler step (n s = 1). The step size e = r/n stcp of the 
molecular dynamics trajectory integration is held approximately constant in all 
data points. (In all runs we targeted acceptance rates between 70% and 90%.) As 
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Figure 1: Integrated autocorrelation time Ti nt versus the trajectory length r tra j in 
CP 9 , (3 = 0.7. r int of the topological charge Q, the magnetic susceptibility xm and 
the energy E is shown. Left side HMC without trivialization, right side THMC 
with one step of Euler integration n s and integration of the flow to tx = 0.47, 
such that the force is minimized. For both HMC and THMC a choice of the 
trajectory length around one is a good compromise. 



already expected from the QCD experience, the optimal value of the trajectory 
length depends significantly on the observable. The energy E decorrelates fastest, 
with a clear minimum at r tra j ~ 0.3, whereas the magnetic susceptibility exhibits 
a very shallow minimum starting from r tra j f=s 1. The topological charge Q can 
profit from even longer trajectories. This is the case for the standard HMC as 
well as the one including the field transformation. Considering the computational 
costs of effectively decorrelated configurations, r tra j in the range between 0.5 and 
1 seems to be efficient. We therefore choose r tra j = 1 for all (3 in our main runs 
as a compromise. Although this might not be the optimal choice, it is a standard 
practice in QCD simulations. 



4.1.2 Integration of the flow 

The second parameter to fix in our setup of the THMC is the value tr to which the 



flow Eq. ( 14 ) is integrated and the accuracy of the integration, which is given by 
the number of steps in the Euler integration. As a criterion, we use the reduction 
of the forces experienced in the molecular dynamics evolution, because perfect 
trivialization would result in forces equal to zero. The relative reduction of the 
force R = Fthmc I 'Fhmc for P = 0.7 and (3 = 0.9 in the CP 9 model is shown in 
Fig. m In both cases a reduction by about 60% can be reached at a value of tx 
around 0.5, much smaller than the tx — 1 for which trivialization is reached with 
the exact flow. The force reduction depends very little on the accuracy of the 
integration: whether 1, 2, 5 or 10 steps of the Euler integrator are used hardly 
matters. As shown in Fig. [3j the optimal value of tx decreases with increasing /3, 
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Figure 2: Reduction of the force in THMC compared to the force in HMC (R = 
Fthmc I Fhmc) depending on the integration length tx of the trivializing flow. 
Data are shown for CP 9 with 1,2,5 and 10 steps of the Euler integration. Left 
side (3 = 0.7, right side (3 = 0.9. At large t T , the curves are ordered as in the 
legend, n s = 1 at the top and n s = 10 at the bottom. The statistical errors are 
too small to be seen. The improvement for n s > 1 is negligible. 

however, the reduction of the force at the minimum is almost constant, at least 
in the range (3 = 0.7, . . . , 1.0 which we have investigated. 

Besides the reduction of forces, the field transformation is also supposed to 
reduce the autocorrelations experienced in the simulation. In Fig. [4] we therefore 
show the integrated autocorrelation time for different flow integration lengths tx 
for CP 9 with n s = 1 and n s = 2. This can also be seen as a check of the force 
criterion depicted in Figs. [2] and [3] As expected, the optimum tx to minimize 
autocorrelations depends on the observable considered. Qualitatively, we observe 
that the force criterion leads to a good choice of tx with respect to reduction 
of the autocorrelation. However, while for larger values of tx the force ratio 
becomes worse, the autocorrelation shows a fairly flat behavior. This test was 
done at fixed molecular dynamics step size and therefore fixed cost per trajectory. 
The main reason for the autocorrelation time rising for large tx is the decreasing 
acceptance rate. Had we kept the acceptance rate constant (meaning increased 
costs for tx larger than the force minimum), the autocorrelation times would have 
decreased further, although not very much. We conclude that the force criterion 
is reasonable for tuning tx and we therefore used it throughout this study. 

4.1.3 Total cost of the simulation 

The reduction in the forces by about a factor two from the field transformation 
allows a larger molecular dynamics step size by about the same factor. In partic- 
ular for larger values of /3, this leads to the same acceptance rate for HMC and 
THMC, see Table [T] for details. In our implementation, an elementary leap-frog 
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t T 



Figure 3: Force reduction R versus the flow integration length tT, only 1 step of 
Euler integration (n s = 1), for several values of in CP 9 . At small tx, the curves 
are ordered as in the legend, = 0.7 at the top and (3 = 1.0 at the bottom. As 
(3 increases, the minimum moves towards smaller tx, however, its depth does not 
change dramatically. 

step of THMC with n s = 1 costs almost three times more than with HMC More 
integration steps will further increase the cost. Together with the increased step 
size, this translates to roughly a factor of 1.5 increased cost per trajectory. In 
the next section, we will find a reduction of the autocorrelation times between 
roughly 1.5 and 1.8, depending on the observable. This means that the total cost 
of simulation for HMC and THMC with n s = 1 are about the same. 

5 Results 

With this setup, we performed extensive runs at correlation lengths between 
£ ~ 2.3 and £ ~ 16.6, using the plain HMC and compare it to the THMC with 
the flow integrated with n s — 1 Euler step. For the latter we use the optimal 
values of the flow parameter t T with respect to reduction of the force. The detailed 
parameters can be found in Table [TJ expectation values of various observables in 
Table [2] The measured autocorrelation times in units of molecular dynamics time 
are listed in Table [3j 

5.1 Critical behavior 

This brings us to our main result, the critical slowing down of the simulations as 
(3 — > oo. For large correlation length £, the autocorrelation times are expected 
to grow as 

r int (A) (X r (20) 
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Table 1: Parameters of our runs in CP 9 with coupling /3 and lattices of size L 2 . n s 
is the number of Euler integration steps, where zero corresponds to the standard 
HMC algorithm. r tra j denotes the integration length of the molecular dynamics 
trajectory, n step its discretization, t T the integration length of the trivializing flow 
and P acc the Metropolis acceptance rate. The last column gives the statistics in 
units of molecular dynamics time. 
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Table 2: Expectation values of our runs in the CP 9 model, n s is the number of 
Euler integration steps, further parameters are found in Table [TJ We give results 
for the correlation length £, the energy, the magnetic susceptibility and the square 
of the topological charge. 
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Figure 4: Integrated autocorrelation time r int versus the integration length tx 
of the leading order trivializing flow, for CP 9 with /3 = 0.7. Left hand side: 
one step of Euler integration n s — 1. Right hand side: n s = 2. As shown in 
Fig. [2j maximum force reduction is obtained by tx = 0.47 (n s = 1) and tr = 0.51 
(n s = 2), respectively. 

with z the dynamical critical exponent. It depends, of course, on the observ- 
able A. This scaling is only expected for asymptotically large £, however, also 
with our limited range we can get an estimate of the severeness of the problem 
and the reduction brought by the field transformation. Since we have periodic 
boundary conditions, topological sectors are expected to form in the continuum 
limit. Because of the ensuing barriers in the free energy, Ref. [5] suggests for the 
topological charge an exponential behavior of the form 

r int (Q 2 ) cxexp(c^) • (21) 

However, as we will see below, the presence of such an exponentially slow mode 
does also have an effect on all observables which do not completely decouple from 
it. 

In Fig. [5] we show the Ti nt of various observables as a function of the correlation 
length for the HMC algorithm without field transformation. As expected, the 
slowing down in the topological charge is much more severe than in the other 
observables which in the interval 4 < £ < 13 show a behavior compatible with 



the power law Eq. (20). Making statistically relevant statements about this is 
difficult, because the acceptance rates are not constant over all our runs. We 
try to compensate for that by considering P aC cT in t, but of course this is only a 
partial correction. Also the data is not expected to follow exactly the leading 
order scaling law; due to the high accuracy of our data next-to-leading orders 
might become visible. Nevertheless, fitting the data in the range of 4 < £ < 9 to 



Eq. (20), we get z = 2.0(1) for the magnetic susceptibility and for the correlation 



length. The errors are statistical only. The energy E exhibits a very flat behavior 
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Table 3: Auto-correlation times corresponding to Table [2} For /3 = 0.95 and 
f3 = 1.0, systematic errors for t(E) account for the uncertainty in estimating the 
contribution of the tail of the auto-correlation function. 



in 2 < £ < 13 with z = 0.12(1). For E and £, the fits have a x 2 /dof between 
1 and 3.3, which is acceptable considering the simple formula and the problems 
discussed above. However, while the behavior for E and £ are compatible with a 
power law up to £ ~ 12.1, the last data point at £ ~ 16.6, and for \m also the 
point at £ ps 12.1, show a clear deviation. We interpret this as a consequence of a 
correlation between these observables and the topological charge and will discuss 
this issue below in detail. 

The square of the topological charge exhibits a much worse scaling behavior 
than the other observables. Fitting Eq. (20) to the data with 4 < £ < 12 , we 
get z ~ 4, however, the agreement is not convincing and the % 2 /dof ~ 2U is 
poor. The exponential function Eq. (21) works much better and delivers a good 
description of the data in the whole region 2 < £ < 17 with c ~ 4.2 and 9 ps U.43. 
It has x 2 /dof rj 0.25, but due to the problems discussed above, this has to be 
taken with care. 



5.2 Effect of the slow modes 

Having detected at least one very slow mode in the simulation raises the question 
to what extent the various observables are affected. The answer will depend 
both on the particular observable and the accuracy required in the simulation. 
We interpret the deviation from the power law scaling behavior of the energy, 
the magnetic susceptibility and the correlation length observed in Fig. [5] at £ ps 
16.6 (and weakly already at £ ps 12.1) to be a consequence of the correlation 
between the slow mode and the observables. As observed in the topological 
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Figure 5: The integrated autocorrelation time of various observables for the HMC 
as a function of the correlation lengths. In this log-log plot, the dashed lines 
indicate the results of power law fits to the data for 4 < £ < 13. The solid line 
represents the the exponential form Eq. (21 ). The error bars are smaller than the 
size of the symbols. 



charge squared, the time constant of this mode rises exponentially and at this 
point its contribution is no longer sufficiently suppressed by the smallness of the 
coupling to the observable and becomes noticeable. 

If we identify, for a moment, the slow mode with the topological charge, 
one can understand the phenomenon with the following: while the simulation 
is trapped in one topological sector, it samples the observable restricted to that 
sector A(Q). If the estimate obtained before moving on to another sector after 
about t(Q) steps is more precise than variance A 2 = vaxq(A{Q)) of A(Q) over 
topological sectors, then the autocorrelation time will have a significant contri- 
bution from the topological modes. What matters thus is the relative size of 
^va,i(A(Q))r- mt (A) / T- mt (Q 2 ) and A. In our data we observe that while the for- 
mer is larger than the latter for most of our data points, this ordering is reversed 
for the points at the largest correlation length. If one is interested in the level of 
accuracy given by A, the simulation has to run over many Ti nt (Q 2 ). 

As discussed in Ref . [2] , the slow modes also pose a problem for the accurate 
determination of the autocorrelation times themselves. By restricting the sum 
in Eq. (19) to some window W, a small in amplitude but potentially long tail is 
neglected. To illustrate this, we show the autocorrelation function of the magnetic 
susceptibility at (5 = 1 for the HMC algorithm in Fig. |6j At the beginning, it 
falls quickly to p(t) w 0.01 but then develops a very long tail, a situation already 
described in Ref. |9j . The tail is compatible with a single exponential with a time 
constant equal to the exponential autocorrelation time extracted from pqiit). 
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This is indicated in the figure by the dashed lines. We can use this information 
for an improved estimate of the autocorrelation time [2]. The usual sum of p(t) in 
only performed up to the point where the single exponential tail starts. The rest 
of the sum is substituted by the integral over the single exponential for which the 
largest observed time constant observed in all observables with the same parity is 
taken. In our situation this is T cxp (Q 2 ). Since the HMC obeys detailed balance, 
this gives a strict upper bound for r int , provided that there are no modes which 
suffer from an even slower evolution. 

Also for the correlation length and at /3 = 0.95 a similar behavior can be 
observed. Even though the coefficients might seem small, the slow modes still 
have a sizable contribution because of the very large time constant. At = 1, 
this tail contributes roughly 30% to Tj nt (£) and 50% to r int (xm); for = 0.95 
the contribution of the tail is roughly 10% and 17%, respectively. The values of 
the autocorrelation time from estimating the contribution of the tail in this way 
lie within the la error of the values given in Table [3] obtained by a summation 
to large values of W. In case of Ti nt (E), the improved estimator is significantly 
higher than the value obtained from the truncated sum. The single exponential 
dominates from t « 150, contributing roughly 10% (50%) at — 0.95 (/3 = 1). 
We take the different values as upper and lower bounds and state the discrepancy 
as systematic error in Table |3j If the true value is close to the upper bound, the 
scaling of T- m t(E) deviates from a power law already at £ ~ 12.1. If we assume 
that the exponential growth in the time constant is not compensated by the 
decrease in the coefficient, this contribution will be even more pronounced when 
(3 is increased further. 

5.3 Performance of the field transformation 

We finally come to the comparison between the HMC and THMC algorithm, and 
we show the reduction in autocorrelation time achieved through the introduction 
of the field transformation. In Fig. [7] we plot the ratio of the autocorrelation of 
our observables for the two algorithms. The correlation length and the magnetic 
susceptibility for which we observe a 40% reduction profit most from the field 
transformation. For the topology roughly a 25% reduction is found, the energy 
is almost unaffected, however, it shows a quite short r int over the whole range 
of data. Note that for £ f» 12.1 and £ ~ 16.6 the reduction of v int {E) has to be 
taken with care due to its systematic uncertainty. All critical exponents of the 
THMC algorithm extracted from the range 4 < £ < 13 agree within uncertainties 
with the ones of HMC. Also the deviation from the scaling law at £ ~ 12.1 and 
£ Pd 16.6 is observed. The exponential behavior of r- lTLt {Q 2 ) is compatible with 
HMC within error bars as well. We can conclude that the field transformation 
does not affect the scaling of these variables in the investigated region. 

As commented above, the improvement factor which we find is close to the 
additional cost of the simulation and therefore the two algorithms perform rather 
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Figure 6: Normalized autocorrelation function for the magnetic susceptibility for 
the HMC algorithm at f3 — 1. The dashed lines correspond to a single exponential 
whose time constant has been extracted from the autocorrelation function of Q 2 . 
The coefficients have been adjusted such that the two lines contain the la region 
around t = 1000. This contribution accounts for roughly half the integrated 
autocorrelation time. 

similarly. These findings do not seem to depend strongly on TV, since we also did 
limited simulations in the CP 20 , with essentially equal results. 

6 Summary 

Whether a modification to an algorithm actually improves its performance is 
often very difficult to predict. It therefore needs numerical simulations to study 
its effects. Here we investigated recently proposed field transformations which can 
lead to a speed-up in HMC simulations. Unfortunately, the result is negative. 
Although we observe a reduction in autocorrelation times, the scaling towards the 
continuum limit is not improved. The reduction in the forces, which can be used 
to increase the step size of the molecular dynamics integration, is compensated 
by the computational overhead of the method. However, this conclusion does not 
have to be universally true for all theories. In QCD with dynamical fermions, 
e.g., the computational cost of the construction would be a minor part of the 
whole cost of the simulation. 

Investigating the pure HMC algorithm serves also as an illustration that ex- 
ponentially slow modes will at some point affect other observables of the theory. 
The deviation at £ 16 from the scaling behavior, which we observe for the ob- 
servables up to a correlation length of 12, is therefore a cautionary tale for QCD 
simulations. Even though the slow mode observed in the topological charge might 
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Figure 7: Cost reduction in terms of autocorrelation time for the different ob- 
servables brought by the field transformation. 

not seem to have any influence on other observables in today's simulations [2], at 
some point, the correlation to the topological charge can also affect the scaling 
behavior towards the continuum limit in these channels. 
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